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1 , INTRODUCTION 


Viscous flows have been calculated for a long time by separating the 
flow into a boundary layer and a potential flow region. The boundary 
layer then Is calculated by using integral methods and the potential flow 
by using inviscid flow equations. This approach does not require large 
computers, and even can be carried out by use of hand calculations. This 
approach, however, has its limitations. For example, in calculating flow 
through nozzles, channels, and diffusers this approach is limited to the 
region near the entrance. As the boundary layer thickens, the flow cannot 
be meaningfully separated into a potential flow region and a boundary layer 
region. Further, this approach cannot describe the corner effects, such 
as the secondary flows, or adequately describe flows in which the body 
forces are nonuniform in the cross-stream plane (such as magnetohydrodynamic 
channel flows). 

In the past two decades, with the advent of high speed computers, 
many numerical methods to compute velocity and temperature fields without 
separating the flow into a boundary layer and a potential flow region have 
been developed. In this paper, a new approach to calculate three-dimensional 
compressible viscous flows is presented. The approach is an extension of 
the author's previous work [1] to three-dimensional flows. The method, 
as presented in this paper, is for "parabolic" flows through rectangular 
ducts, and its possible extensions to other types of three-dimensional 
flows is under investigation. "Parabolic" or "boundary layer" flows are 
characterized by the existence of a predominant flow direction along which 



downstream conditions have a negligible Influence on the upstream conditions. 
This assumption allows use of a marching Integration procedure. A lucid 
discussion of the parabolic flow assumption Is given by Caretto et al [2], 
and Patankar and Spalding [3]. 

Methods developed so far to compute three-dimensional viscous flows 
are, among others, those given In [2-10]. These methods differ In their 
finite difference approximations, and whether the flow Is assumed to be 
parabolic or not. From among these methods, those proposed for calculation 
of confined flows differ further In the procedure used to calculate the 
pressure field. These methods are discussed here briefly. . Included in the 
discussion Is the method of Harlow and Welsh [11]. Although this method is 
for two-dimensional confined transient viscous flows, It has some features 
In common with the methods developed later for computation of three-dimen- 
sional parabolic flows. The method of Chorin [4] calculates incompressible 
elliptic flows. Finite difference equations are obtained by using the leap- 
frog method for the time and convective terms, and Duffort-Frankel for 
the diffusion terms. The pressure field Is calculated by coupling the 
incompressible continuity equation to the momentum equations via an artifi- 
cial variable density and an artificial equation of state. Miller [5], on 
the basis of the work of Chorin, computed steady Incompressible three- 
dimensional elliptic flows through rectangular ducts. In the method of 
Harlow and Welsh [11] the finite difference equations are obtained by 
using forward in time and center in space differencing, The pressure 
field Is calculated from an elliptic equation obtained from the con- 
tinuity equation and the divergence of the momentum equations. The 
method uses a "staggered" grid in which the velocities are defined at 
the cell boundaries and pressure at the cell centers. A similar type of 
staggered grid is used in the cross-stream plane, later, in [2,3,6, 8] for 
the computation of three-dimensional flows. In these methods the cross- 


stream velocities are specified at the cell boundaries and the axial veloc- 
ity, along with all the other variables, at the cell centers. Methods 
[2, 3, 6, 7] calculate three-dimensional confined parabolic flows. In these 
methods the pressure field In the cross-stream plane is calculated from 
an elliptic pressure equation obtained from the continuity equation and the 
cross-stream momentum equations. The axial, along the duct axis, pressure 
gradient is calculated to Insure conservation of mass along the duct axis. 

In these methods Caretto et al. [2], and Patankar and Spalding [3] obtain 
finite difference equations by using the control volume approach; Emery et al. 
[6] use Dufort-Frankel for the axial velocity momentum equation and implicit 
differencing for the cross-stream velocity momentum equations; Roberts and 
Forrester [7] use one-sided upstream differencing in the flow direction and 
second order differencing in the other two directions. Partap and Spalding 
[8] extend the method of Patankar and Spalding to "partially parabolic" 
flows. Methods of Nash [9] and Wang [10] are for external (pressure field 
specified) three-dimensional parabolic flows; Nash uses forward explicit 
differencing, and Wang uses Crank-Nicolson type differencing scheme. 

All the methods described in the preceding paragraph use Eulerian 
system. Flow field is computed by calculating the velocity components 
u v , u u , and u along a set of spatial grid points fixed in advance. In 
contrast, in the present approach, the flow field is computed by calculating 
the streamwise (not axial) velocity, u, along a set of chosen streamlines. 
Streamlines are identified by a pair of indices, (i,,j), such that the 
coordinates of the streamline (i,j) in the cross-stream plane are and 
Z.. Coordinates of the streamlines in the cross-stream plane are calculated 

J 

as part of the flow computations. The dependent variables in the present 

approach are u. ., the streamwise velocities, and ( Y . ,2 . ) , the coordinates 
• *J * J 



of the streamlines. To calculate the streamwise velocities streamtubes are 
constructed around the individual streamlines such that the sum total of 
all the streamtubes fills the duct. One such streamtube, constructed 
around the streamline (1,j), is shown in Fig. 1.1. Finite difference equa- 
tions for u^j and h^j, velocity and enthalpy along the streamline (i,j), 
are now obtained by applying Euler's momentum theorem and the first law 
of thermodynamics directly to the flow in the streamtube between x and x + Ax. 
After the streamwise velocities and densities have been calculated at x + Ax, 
the streamline coordinates at x + Ax are calculated next by using mass con- 
servation. From (Y-j.Zj) and (Y^zj 1 ), coordinates of the streamline (i,j) at x 
and x+ax, one can calculate the slope of the streamline (1,4), and then using 
this slope decompose u. into its components u 4 x u/ and uf <. 

1 1 , j * >j i , j i * j 

Strictly speaking, the present approach is not Lagrangian, as one does 
not follow the individual fluid particles, but rather follows the flow 
through individual streamtubes. This approach, however, seems closer to 
the Lagrangian than the Eulerian flow description. The advantage of the 
present approach is that one needs to solve only one set of finite difference 
equations— for the streamwise velocity u— rather than three sets of equations 
for u x , Uy, and The disadvantage is that extra computations are required 
to calculate streamline coordinates. 
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2. FINITE DIFFERENCE EQUATIONS 


2.1 Background 

Consider the flow through a rectangular duct with Its axis along the 
x-dlrectlon. Flow In the duct Is partitioned Into a finite number of stream- 
tubes. To define the various streamtubes lines Y I are 

drawn parallel to the z-axis (left duct wall), and lines Z2, Z 3 ,,.,, Zj,..., 

Zj parallel to the y-axls (lower duct wall), as shown in Fig. 2.1. Lines Yj 
and Zj (not shown In Fig. 2.1) are drawn through the middle of the lower 
wall and the middle of the left wall, respectively. The lower and left 
Walls are also referred to as Z\ and Yj lines, Line Y^ is drawn in the 
middle of Yj and Y 1+1 , lines Y^* Zj +i , and Zj.$ are drawn similarly. 

Various streamtubes, called streams for short, are now defined as follows. 

Stream (2,2): flow bounded by the left wall, the lower wall, Z 2 + j, 
and Y 2 +j. 

Stream (1,2) : 1*3,1: flow bounded by the lower wall, Z 2 +J, Y^.^, 
and Y 1+i . 

Stream (2,j): flow bounded by the left wall, Y 2 +*. Zj. if 

and Z j+i . 

Stream (i,j): 1*3,1; j*3,J: flow bounded by Y^.|, Yj+$, Zj_|, 
and Zj+j. 

Streamline (i,j) is defined as the streamline passing through the intersection 
of Y.j and Zj. 
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the various streams. 


Other nomenclature 1$ as follows: 


fjji mass flow rate through stream (i,j). 

The mass flow rate yjj Is partitioned further Into four parts y*j, 

and mass flow rates through the upper-right, upper-left, lower- 
left, and lower-right quadrants, respectively (see Fig. 2.1). 
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V 

j: cross-sectional area of the stream (i,j). 


To apply conservation principles control volumes are constructed for each 
stream using a distance Ax along x. Control volume for the stream (1,d) is 
shown in Fig. 2.2. 


A? 1 ; Af J : surface areas of the control volume normal to the y and z 
j i 

axis respectively. 



T 1,j*’ N,J* p i,j’ c i,j ; velocity, temperature, enthalpy, density 

and constant pressure specific heat along the streamline 
(1»j) at x. 


ut ., etc.: variables along the streamline (i,j) at x + Ax, 

* 

u 1tj : average of j and ujj. 

u? <1 u^f u? components of u . . along the x, y, and z directions. 

iij 1 t J * * tJ 




Yji Z j : coordinates of the streamline (1,4) In the cross-section 
plane at x. 

Y*; zj: coordinates of the streamline (1,j) at x + Ax. 

S U 2 *“i.j 

s u 5 * “u 

^l^.d-i 5 e ^ ect ^ v * viscosities, Including the 

turbulent contributions, for the momentum transfer across the 
surfaces of the control volume, as shown In Fig. 2.2, evaluated 
using the Information at x. 

r i 

w i+i j* etc ‘ ; e ff* c ^ve viscosities evaluated using the Information 
at x + Ax. 

K i+i j’ * tc,: e ^ ee *i ve thermal conductivities, defined similarly as 
the effective viscosities. 
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4. 

Corresponding V's and H’s with superscript + are defined in terms of p , 
k + , Y + , and Z + . 

In estimating momentum and energy fluxes through a streamtube u and h 
are assumed, in general, to vary linearly between the neighboring stream- 
lines. Between the walls, however, and the adjacent streamlines (streamlines 
with i or j equal to 2 in Fig. 2.1), a nonlinear variation of u and h is allowed. 



Thus, the mass, momentum, and energy flux through the streams adjacent to 
the walls are expressed using coefficients y, a, 0, and 6 defined as follows 
Mae Appendix A for further discussion): 


h f v 2 


a ‘I f 


pudy dz 5 y* pu 2>2 Y 2 Z 2 


0*0 

Z 2 


n 

* o J o 

Zo 


pu u dy dz = 


_ z 


I 

•' 0*0 


pu i \i c dy dz 5 


“1 V 2,2 u 2,2 


6 1 V 2,2 * u 2,2 


f *2 f Y 2 

I j Puhdydz 5 0*f 2 ; 2 h 2f2 

» n J n 


The following relations are for the lower quadrants of the streams along 
the lower duct wall. Formulae for the lower-left quadrants are obtained 
from these relations by substituting for 2 , 


aY 


Y 1 “ Y M - _ 3u 1,2 + u i-l ,2 

~ — 2" * U 3 -1 » 




P 3iu^ 2 + iuAi 2 
and iu z = — -*—4— - - 


and for the lower-right quadrants by substituting y^” 2 for ^ 2 


AY 


Y i-H 2 “ Y i f , 3lJ i ,2 +^1+1,2 


3h i ,2 + h 1+1 ,2 


an« . 3>U <.2 - 2 . 
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*1.2 * AY 


f Z2 

f I p u dz 

•'o 

f L 2 

AY I puudz 1 

J 0 

C h 2 

AY I pu $u* dz s 
«'0 

c lz 

AY I p u h dz s 

J n 


Y* p U AY z 2 , 






0) z y 1i2 K . 


Coefficients yj» ct^, frij, and 3j» that take Into account the non-linear 
variations normal to the left wall, are defined similarly. 


2.2 Finite Difference Equations for u^ j 


Applying Euler's momentum theorem to stream (i,j) between x and 
x + Ax one obtains. 


Momentum flux out 


at x + Ax 


- Momentum flux ini 

[at x 


* Pressure forces + Viscous forces + F. . (Body forces) 

* 

(2-1) 

For the streams (i,j) with i and j>2, various terms in Eq, (2-1) are 
estimated as follows. 

Momentum flux (in or out) through the stream is calculated by summing 
the momentum fluxes through each of the four quadrants (see Fig. 2.1), 

Momentum flux through any one quadrant is estimated as equal to the mass flow 


rate through that quadrant times the average velocity, which is obtained by 
double Taylor series expansion of u^ j around i,j, for that quadrant. One 
thus gets; 


Momentum flux 



- [u + *u]. 


out 


- Momentum flux in 


I at x + Ax 

+ u 1+l .j + “U+l ) 




at x 


( z Vj + U u+i 

( 2u . + . + u t . . 
\ i.j V.J-1 



( 2 - 2 ) 


± j. 

Symbol [u + u] denotes the quantities in the preceding bracket [ ] with u 
replaced by u. This symbol is used similarly later on. 

In calculating pressure and viscous forces, it is assumed that the 
streamline inclinations with respect to the duct axis are small; that is, 
t/ and u z are small compared with u x . Thus, in calculating pressure force 
A, x j , the cross-sectional area of the streamtube (i,j) projected onto the 
Y-Z plane, is used. In calculating viscous forces, the velocity gradient 
is based on An x (see Fig. 1.1), the distance between the neighboring stream- 
lines projected onto the Y-Z plane. Pressure force on stream (i,j) is 
estimated as 

- -DPA^j (2-3) 

where DP is the pressure difference between x and x + Ax and is assumed to 
be the same for all the streams. Viscous force on a surface between x and 
x + Ax is approximated by the average of the viscous forces evaluated by using 


the effective viscosities and the velocity gradients at x and x + Ax. One 
thus obtains for the viscous forces on the stream (1,j) (see Appendix A), 


Viscous Forces - \ [v^.j <V),J ’ <,j> + ^, + J+l <u 1 + ,J*r u U> 

* v 1 I .V“U' u l + -l.J ) ' v u (u l.j ' u 1 + .j-l ) ] 


+ \ [V + ,u+-V,u], 


(2-4) 


For streams along the walls, momentum fluxes are expressed using the 
non-linear correction factors introduced in section 2.1. Viscous forces 
exerted by the lower and the left wall on the adjacent streams are, respectively* 
set equal to 

V 1 J ,2 u i,2 and ^2 , j u j ,2 

1 J 

Calculation of Vj^ 2 and V 2 j is discussed in Appendix A, 

After the various estimates in Eq. (2-1) have been substituted, all 
the terms Involving u^ + j are collected on the left-hand side* and the resulting 
equation Is organized in the form: 


- A, ^ u . + , i , - C / * u. + , . + B . . u.+ . - A. J . u, + . . - C, J . u * , 

1*3 1 + 1 *J 1 ij 1 > J 1*0 1 t3+l ^ t J 


D i,o 


(2-5) 


Coefficients A, B, C, and D for the various streams are given next. 

NOTE : In equations that follow all the Vs, except VjT'g, are * tlnles 

the f s defined in section 2.1, and the V's are i times the 
Vs defined earlier. 



Equations (2-6) : 


Coefficients for stream (2,2) 


A 1 * - 4f ++ 

* 2,2 * 2,2 


°2 '*' 2, 2 + V 3 J , + 2 


C 2,2 * 0 


2.2 ■ + 3a 2 <2 * «f V.* + 3a 2 + *3,2 + »2* 3 ♦ V& ♦ V 


J + 

2,2 


-f 0 ++ , - a? + V 


J+ 


2,2 2,2 * a 2 r 2,2 r v 2,3 


C 2 ,2 " 0 


' 2,2 


“ DPA 2,2 + F 2,2 


U 3,2 I ^2+2 + a 2 ^2^2 * V 3^,2 


+ U 


2,2 


v. 


2 f 2,2 + 3a 2 ^ 2,2 + + 3ot 2 V 2 + ,2 


V 3,2 V 2 J ,3 " V 2,2 " V 2,2 


U 2,3 + w 2 ¥ 2, + 2 + V 2 J ,3 ) 
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Coefficients for stream (1,2), 1>2, 


A /,2 " "'*' i + , + 2 " a f **,2 + V 1 +l ,2 


C 1 X ,2 * <2 " “f *i72 + V i, + 2 


B i ,2 " 2H, i°, + 2 + 3a 1 Z V UZ + V i+ + 1 ,2 + V M + V 1 I , + 2 + V i°, + 2 


A 1 J ,2 * ~*U2 + V 1 J , + 3 


Cfe ■ 0 
* ,*• 


°1,2 - - DPA Z + F 1,2 


Vl ,2 ( V 1 + , + 2 + “t ' l 'l + ,2 + C .2 ) 


+ “1-1.2 (’Ca * “* V.2 * *f 1 .2 


+ u i ,2 (< + 2 + 3 “1 2 C + Cl ,2 + < + 3 + V l‘, + 2 + < + 2 


+ U 1 + C) 


Coefficients for stream (2,j), j>2. 


A 1 . -Y+ 0 + v I+ 
A 2,j ’ V 2,J * V 3 J 


C 2.j " 0 


•W ■ + Hf ’iS + »3 , .J * »&♦! * + »& 


■ .ur"*"* . a y w 

A 2.j *2,j a j f 2,J + V 2,j+1 


c ^ » -«y u/ _ u» + v 

L 2,j a j t 2,J t 2,J + V 2»j 


D, 4 » -OP A* + F 0 4 




2J 


u 3,d(' 1 '^. t j + V 3,j) 


+ U 2,j (»2J + 3 “j *2J ‘ »J.j ■ V 2.j + 1 ’ V 2,j - V 2 J 
+ U 2J + l(< + j +V 2 J ,J + l) 


U 2,j-l(“jV;j + <‘4 +V 2 J .j) 
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Coefficients for stream ( 1 ,J), 1 and J> 2 , 


. I _ v +0 . y I+ 
A 1 J “ f 1 .J + VlJ 


. I . » -0 . y 1 + 
•i.d -v t.j + v i.j 


B 1 ,i * *Ui * V 1+' ,J + W 1 J .V' + V u + V ^3 


id . -V 0+ + V J+ , 
'l.J *1.J V 1 .4+1 


. 0 * w 0" + y d* 

'l.J ¥ i,J v i.J 


0 . , * -DP A/ h * F, , 
l,j l.J !»J 

+ u 1,j( T l“ ‘ V H1.J ■ V U+1 ' V U ' V < J .^ 

♦ (*u + + ( f u + v '^ 
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2.3 Finite Difference Equations for h^ j 

Applying first law of thermodynamics to the stream (1,j) between x 
and x + Ax, one obtains, 

Energy flux out - Energy flux in * 

at x + Ax at x 

Heat added by conduction + Work done by the viscous forces 
+ Qj j( Internal heat generation) + W^ j(Work done by the body forces) 

(2-7) 

For the streams (i,j) with 1 and j>2, various terms In Eq. (2-7) are 
estimated as follows: 

The energy flux terms are estimated In the same way as the momentum flux 
terms and can be obtained by substituting (h+s) for u In Eq. (2-2). 

Heat conducted across the top surface of the control volume (see Fig. 

2.2) Is estimated as 

“ ?[ H fJ+l <h 1 + ,o+l ■ h i + ,j>] + I [hW-’-M] (2-8) 


Work done by the viscous forces on the top surface is taken to be 



(2-9) 


Work done by the viscous forces and the heat conducted across the other 
three surfaces of the control volume are estimated similarly. 
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After the various estimates In Eq. (2-7) have been substituted, all 
the terms Involving h^ + j are collected on the left-hand side and the resulting 
equation Is organized In the form: 








( 2 - 10 ) 


Coefficients A, B, C, and D for the various streams are given next, 


NOTE : In equations that follow all the Vs, except H' 2 *" 2 , are } times 
the 'tf's defined In section 2.1, and the H's are i times the 
H's defined earlier. 


Equations (2-11) 


Coefficients for stream (2,2) 
* 2^2 " m * Z *2 " 4 * 2,2 + H 3, + 2 


C 2,2 * 0 


2.2 ■ 2f 2?2 + K ha * *f *,-,2 ♦ * H ,»♦ ♦ H 3 * ♦ H & * M* 


*2.2 ■ <2 ha * <3 


h 3,2 

(h*a * h z h*a + <2 

) 




h 2,2 1 

( 2 ha * 3 h y h'a ♦ B f 

T 2'.2 * 

30| 

w + “ 
y 2,2 

" H 3 ! ,2 ’ 

^2,3 ( 

h'a ♦ s 2 <2 ♦ h 2 j >3 

) 




S 3,2 ( 

> 2*2 + 4 ha + h : a 

) 




S 2.2 ( 

2 h*a + 33 2 f 2, + 2 + «f 

h'a + 

3«| 

m+- 

*2,2 

’ V 3,2 “ 

s 2,3 1 

( y 2,2 + s 2 ^2 ,2 + V 2,3 

) + h i 

,2 H 2 

‘ 2 - 

h 2,l H 2 J ,2 

I ’1.2 H 2*2 + h 2,l H 2 J , + 2 + s 3,2 ( 

4,2 ’ 

*1 

y/ + “ 
*2,2 

* * ) 


+ s 2 + ,2 (» 2?2 + 33 2 T 2.2 + 6 1 h'a + 3 i 2 * <2 + V 3.2 f V 2, + 3 + W 2. + 2 + V 2 , + 2 ) 
+ S 2 + ,3 (< + 2 - «? <2 ♦ »&) 
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Cotfflcltnts for strtam (1,2), 1 > 2. 


t I m »++ _ a 2 w+- * u ^ 

* 1,2 M ,2 S 1 M ,2 + H 1 * l ,2 


C l ‘,2 ■ <2 ■ 9 1 * l '.*2 + H /*2 


*1.2 ' < + 2 + 38 1 *“'2 + H H1,2 + < + 3 ♦ <2 * <2 

A J , ||, 0 + . yj + 

A i.2 " f 1,2 + H 1,3 


C 1 ,2 " 0 


° 1 ,2 * Q 1,2 + W 1,2 

* h 1+l,2 {*UZ * ®1 V 1*2 + H /+l ,2 ) + h i-l ,2 (*UZ + e f + H /,J 

* "|.3«Z + 36 1 2 *1% - H 4l,2 - <3 ■ <2 - <2) 

♦"l^Wj 

+ *1+1 .2 ( .2 + S 1 T , ,2 + V 1+1,2) + s 1-l,2 (' l ',, + 2 + 4 f T 1 ,2 * ' / t I ,2 ) 

+ * 1.2 ( <2 + 3S 1 <’2 ‘ * 1 + 1,2 - V 1 J ,3 ’ ‘‘,2 - < 2 ) 

+ S 1 ,3 ( * V 1°,3 / + h 1,1 H 1 J ,2 + h i ,1 h 1 J , + 2 

+ *1+1.2 \ -< *'l .2 ■ S 1 Z .2 + V i+1,2/ + *1-1 ,2 (■' l 'l, + 2 " S 1 .2 + V 1*2 . 

+ S 1 *2 24, 1°,2 + ^i 0 f 2 + V i+1 ,2 + V i\3 * V /, + 2 + V U2/ 
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Cotfflcltnts for str«*m (2,J), j>2. 




* ul+ 

•*2.J + H 3,J 


I 


C 2,J * 0 


B 2J ■ «2?J * »? ♦ H** ♦ ♦ H «♦ ♦ K" 

A 2J ‘ ’’Cj ’ B J + H 2,Vl 


C 2,j ' * B j f 2~.) ' f 2 + ,'j + H 2 J , + j 


'l.j * 9 2,J + W 2,j 

+ h. 


+ s. 


V 2J + 

H 3.J 

) 




zv z!j. 

+ 3B J 

ur -0 

*2.J 

■ H 3,j ■ 

u J 

H 2,j+1 ’ 

H 2,J ' 

('& 

+ 9 y 


+ H 2 J ,J+1 ) 

^ + h 2.J-l 

I 8 ! *2', 

».V 

V 1 ' 
V 3,j, 

) 




Zf 2J 


m-0 

V 2,J 

„ w I 
V 3,j 

V J 

V 2,j+1 

V I 

V 2J 

( '£j 

+ «? 
j 


+ V 2 J ,j+l ) 

1 + S 2.M 

( s 3 *S 


fj 


+ "iXj + "X.W.J 


+ *3.i K° * ^ 


- S 2,j«j + 3«J *2J + »3. + J + »2J*I + V 2.J + V 2. j ) 


/ ++ y -+ 

+ $2.1+1 ' -Vo 4 - Vo 4 


J+ 


S 2,j+1 \^2,j “ 6 j + V 2,j + 1 / + s 2,j-l\" 6 j *2J- *2,j + V 2J/ 


f y 


j+ 1 
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Cotfflcltnts for strum (1,J), 1 and J>2. 


A /.j ■ <i ♦ H tVi j 
c u ■ -*u * h m 


D B «)00 * |ij,+ , x U 1 + i 

M.J H 1+U H U+1 + M 1,J + H 1,j 

A J m _U/ X U 

*1,j 


c u ■ <i + <i 


’1,J 




+ h l+1 .J + H 1+l ,j ) + h 1-l J ( ’’u * H u) 

+ h 1,j (*“ - m h - h u-<j) 

+ s 1+l,j (»U + v 1+l ,J ) + *1-1 »j ( + v /,j ) 

+ s i.j(<V v iV, >j -< j+ , - v u • v i J ,j ) 




+ s 1+l .J l ' y 1 + .°J + V 1+i .j ) + s i + -l ,j ( + V 1J 


<i ( *“ ♦ c, ,j * 


+ *id + i(- T u + <Vi) + s i + ,j-iK‘j + <j) 



2.4 Decomposition of u and the Secondary Flows 


Computation of secondary flows, defined as the velocity components normal 
to the main flow, Is discussed In this section. After Y* and Zj + have been 
calculated, the streamwlse velocity u^ j, the average of u^j and u^ j, 
then can be decomposed Into ujj, ujj* and u 2 j, by using the formulae, 


ir 


“i.j 


Ar 


and 


iL. & 
Ar 


( 2 - 12 ) 


where Ax Is the Integration step. 


AY-Y 1 + -Y i , AZ * Zj + - Zj » and Ar 2 * AX 2 + AY 2 + AZ 2 


The calculation of Y*and Zj*" as discussed In Appendix A Is based on 
the conservation of mass alone. Velocity components u y and u 2 , obtained 
from Yj^and Zj + thus calculated, are the secondary flows associated with 
the flow development in the ducts— may the flow be laminar or turbu- 
lent, In turbulent flows through rectangular ducts, besides the primary 
shear stresses which are In the streamwise direction, there are also shear 
stresses in the plane normal to the streamwise direction (see, for example. 

Launder and Spalding [12]). Force due to these stresses will produce bending 

+ + 

of the streamlines in the z-x and z-y planes, and thus displace Y^ and Zj 

beyond that displaced by the conservation of mass. For example, force per 

unit area Ax , due to the cross-stream shear stress x distribution, 
yz yz 

will bend the streamline (i,j) according to the formula, 



Ax 


JK. 


Y i + * 


V, 


M 


(2-13) 
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where R^ j is the radius of curvature in the z-x plane. With R* j (approxi- 
mated as a straight line for small Ax) known one can calculate displacement 
of Zj + (and similarly of Y^) caused by the cross-stream turbulent shear 
stresses. Contributions to u x and u*, as calculated from Eq, (2-12), by 
this additional displacement of the streamlines In the cross-stream plane 
represents, as viewed from the present computational approach, the secondary 
flows associated with turbulent flow through noncircular ducts. It should 
be noted that the turbulent shear stresses In the cross-stream plane do 
not enter into the finite difference equations for the streamwise velocity 
presented In section 2.2. These equations are valid for laminar and tur- 
bulent flows as long as the streamwise shear stresses, appearing In the 
equations via V's, are modeled appropriately. 


4 



3. SOLUTION TECHNIQUE 


A methodology to solve the momentum and energy finite difference equa- 
tions given In section 2 Is presented in this section. This is the method- 
ology used for the sample flow computations presented In section 4. Equa- 
tions (2-5) and (2-10) are of th» general form: 


■ A U X 1+1J ' C 1 I ,J X 1 + -l,j + B U X i U 


A ^ X + 

A i.j+n,j+i 


C J v + 

s.ru-i 


= D 




(3-1) 


Equations (3-1) constitute a set of coupled nonlinear algebraic equations. 
The solution technique discussed is for flow through a rectangular duct, with 
the duct cross-section specified along the axis. In this case, the pressure 
drop, DP,, that appears in ^ for the velocity equations, is not known in 
advance and is to be calculated as a part of the computations. Equations 
are solved iteratively by using two iteration loops, one placed inside the 
other. The outer loop iterates on the unknown DP, and the inner loop 
solves the nonlinear algebraic equations, given a value for DP, iteratively. 

The following steps are for one integration step from x to x + Ax. 


1. Assign a value to DP. At the beginning of the outer iteration 
loop, the value assigned is based on the upstream value of DP. 

For the subsequent iterations, DP is assigned a value as discussed 
in Step 7. 

2. Calculate coefficients A, B, C, and D for the velocity equations, 
Eqs. (2-5), evaluating V's by using the information at x for the 
first iteration, and the information at x + Ax obtained in the pre- 
vious iteration for the subsequent iterations. 
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3. Once A, B, C, and 0 are known Eqs. (2-5) constitute a pentadiagonal 
set of linear algebraic equations. These equations are solved by 
using alternate direction Implicit method. 

4. Repeat steps 2 and 3 for the enthalpy equations, Eqs. (2-10). 

5. With the newly calculated values of dynamic and thermodynamic 
variables at x + Ax, calculate Y^ + and Zj + by using the formulae 
given In Appendix A. 

6. Check newly calculated u* j and h^ for convergence. If the con- 
vergence check Is met, proceed to the next step; otherwise repeat 
the steps starting from 2. 

7. Check Yj + and Zj + , calculated duct half width and half height for 
the value of DP assigned In Step 1, against the duct dimensions 
specified at x + Ax. If the check Is within the desired accuracy, 
the computation of flow from x to x + Ax Is complete. If not, DP 
Is adjusted, on the basis of the difference between the specified 
and the calculated duct dimensions, and the calculations started 
again at Step 1, with the new value of DP. The formulae used to 
adjust DP for the succeeding iteration Is the one given in [1]. 
This formula, based on the difference between the calculated and 
specified duct cross-sections, is applicable to both two and three 
dimensional flow computations. 

Spacing between the streamlines in the cross-stream plane changes from 
x to x + Ax. As a result, at the beginning of each new integration step, 
strear^tubes ire defined anew in accordance with the definitions given in 
section 2.1. 


4. SAMPLE FLOW COMPUTATIONS 


The approach, presented In this paper, to compute three-dimensional flows 
by calculating velocity and enthalpy along the streamlines, was checked by 
making three sample flow computations. All the computations are for laminar 
flow through rectangular ducts. Computation of laminar flows Is free from 
the uncertainties associated with the modeling of turbulent shear stresses 
In three-dimensional flows and, thus, provides a clearer assessment of the 
numerical method Itself. 

All the computations are for the lower-left quarter of the duct. Each 
sample flow computation was carried out using one hundred streamlines; thi * 
corresponds to, Including the grid points along the walls, a 11x11 grid 
In the y-z plane. Streamtubes around these streamlines were redefined at 
the beginning of each new Integration step, as discussed In Section 3. 

In the present computations, however, the same set of streamlines was used for 
the entire length of the duct. The integration step along the duct. Ax, was 
selected for each sample flow computation so that the locations of the com- 
puter output will match the locations of the reported experimental data; 
the integration step varied from one-tenth to one-twenty-fifth of the equi- 
valent duct diameter. Velocity distribution at the inlet In all the sample 
flows was taken to be uniform across the duct cross-section. The following 
nomenclature is used in the figures in this section: 

— • — • — : Dots along the velocity distribution curves represent 

the streamline locations along which the flow was computed. 
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U: Axial velocity. 

U Q ; U at the Inlet. 

U c : Centerline velocity. 

W;H: Duct half width and half height, respectively. 

D: Equivalent diameter (» 4WH/(W+H)). 

X: Distance along the duct measured from the inlet. 

Y and Z: Distances normal to x measured from the lower-left duct corner. 
Re: Reynolds number (■ UqD/v, where v Is the kinematic viscosity). 

Pr: Prandtl number (■ 0.72). 

Gr: Graetz number (* PrRe/x). 


Num: 


Logarithmic mean Nusselt number 



T - T 

center wall 



T bulk : tem P eratl < re * 


In Figs. 4.1a, 4.1b, 4.1c, and 4. Id the results for flow through a 
square-sectioned duct are shown. Calculations are compared with the experi- 
mental data of Goldstein and Kreid [13] (Figs. 4.1a, 4.1b, and 4.1c) and 
Beavers et al . [14] (Fig. 4. Id), Agreement between the theoretical and experi- 
mental results is good. In Figs. 4.2a, 4,2b, 4.2c, and 4. 2d the results for 
flow through a rectangular duct of 2:1 aspect ratio are shown. Calculations 
are compared with the experimental data of Sparrow et al . [15]. Once again 
agreement between the theoretical and experimental results is satisfactory. 

In the two flow computations just described only the finite difference 
equations for velocity were solved. Density was considered constant and, 
thus, the energy equation was bypassed. To check the finite difference equa- 
tions for enthalpy a sample flow with heat transfer to the walls through a 
rectangular duct of 2:1 aspect ratio was computed, The inlet temperature and 
velocity distributions were taken to be uniform across the duct cross-section. 




Fig. 4.1a. Development of the velocity profile In a square-sectioned duct— diagonal plane. 







Fig. 4.1b. Development of the velocity profile in a square-sectioned duct— central plane. 


2.2 
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Fig. 4.1c. Center line velocity development in a square-sectioned duct. 







Fig. 4.2b. Development of velocity profile in rectangular duct of 2:1 aspect ratio— across short side 



Fig. 4.2c. Center line velocity development in rectangular duct of 2:1 aspect ratio. 




Fig. 4. 2d. Pressure drop in rectangular duct of 2:1 aspect ratio. 



The walls were maintained at a constant temperature of 100 degrees Celsius 
below the Inlet flow temperature. The velocity and enthalpy finite difference 
equations were solved together as discussed In Section 3. Calculations were 
carried out using constant specific heat, perfect gas behavior, and the 
Prandtl number equal to 0.72. Results of the computations are shown in 
Fig. 4.3, and are compared with the calculations of Wlbulwas as reported In 
[16] Table 52, p. 220. 

To Illustrate th.? velocity field In the cross-stream plane. In Fig. 4.4 
are plotted the streamline locations In the lower-left quarter of the duct 
for air flow through a 3/2 meter by 3/4 meter rectangular duct (second sample 
flow computation) at two locations along the duct. The beginnings and the 
ends of the arrows represent, respectively, the streamline locations at 
1 meter and 5 meters from the duct Inlet. Bulk velocity through the duct 
Is 0.15 cms/sec. 
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Fig. 4.4. Motion of the streamlines in the cross-stream plane- 
figure is for the lower-left quarter of the duct. 



PRECEDING PAGE BLANK NOT FiLMcD 
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APPENDIX A 

In this appendix formulae which supplement the finite difference equations 
given In Section 2 are presented. The lower-left corner of the duct Is 
taken as the origin of the y-z plane. 

A.l Mass Flow Rates 

Let R 1.j ■ 

Formulae for Stream (2,2): 

^ 2,2 " Y 1 R 2,2 Y 2 Z 2 

f ++ „ 2R 2,2 * R 3,2 + R 2,3 | Y 3 ' Y 2 | | Z 3 ~ Z 2 | 




Formulae for Stream (i,2), i>2: 
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Formulae for Stream (2,j), j>2: 



Formulae for Stream (i,j)» i and J > 2 : 



A. 2 Surface and Cross-Sectional Areas 




A. 3 Streamline Coordinates Y*, Zj + 

After the velocities and densities have been calculated at x+Ax, 
the following procedure is used to calculate the streamline coordinates at 
x + Ax. 

Let R u ■ p i\j U 1 + .J • and 



For a given value of j, by summing the mass flow rates through the stream 
tubes with i * 2 to i * I, one obtains for j * 2, 
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and for j > 2, 


< z j- Z j-1 )Y I + 


16 ¥ 


u_ 


16 ¥ 


-+ 


Ll±L 


Y j (3R 2 + ,j + R 2 + ,j-l* Y j y -1 * 3R 2 + ,j-l + R 2 + ,j) 


+ 16 


1-3 


(A-2) 


By adding Eq. (A-l) to Eqs. (A-2) for j - 3 to j»J, one gets 

Z X * Si - (A - 3 

where S: represents the sum of the right-hand sides. By using the duct 
aspect ratio, specified at x + Ax, Eq. (A-3) can be solved for Yj + . With 

X 4* 4* 

Yj known, Eq. (A-l) and Eqs. (A-2) can now be solved for Z 2 , Z 3 , . .., Zj . 

XX 4. 

The coordinates Y 2 , Y 3 , .... Yj are calculated similarly. 


A. 4 Mass, Momentum, and Energy Correction Factors 

For laminar flow, the streamlines adjacent to the duct walls can be 
chosen near enough to the walls so that the velocity and enthalpy distri- 
butions between the walls and the adjacent streamlines can be approximated 
by linear functions. In this case 2 and V 2 y introduced in Section 2.2, 
are given by the general expressions for V's given in Section 2.1, and the 
correction factors y, a, and 6 take on the following values: 
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These are the values for y, ot, and <5 used In the sample flow computations 
given In Section 4. The finite difference equations for the enthalpy 
were solved by using the wall temperature as the reference temperature; 
for this case the 0's become equal to the corresponding a’s. 

For turbulent flow the correction factors y, a, 6, and 0, and V.^ 9 
and ^2 j will depend on the turbulence model chosen for the near wall region 
In three-dimensional flows. Their calculation can be carried out by using 
three-dimensional extension of the procedure given in [1] for the case of 
two-dimensional flows. 

A. 5 Viscous Forces on Stream (i,j) 

The viscous forces exerted on the stream (i,j) by the neighboring streams 
are calculated as follws. The force exerted on the stream (i,j) by the 
stream (1,j+1). call it F(top), is given by 

F(top) = Ayf| • (A-4) 

where A is the interface area between the streams (i,j) and (i,j+l), and 
n is the effective viscosity at this interface. By substituting the appro- 
priate symbols for A and y defined in Section 2.1 (see Fig. 2.2), and by 
approximating the velocity gradient by 



one obtains from Eq. (A-4) 


«**> * <«) 

By using the definition of Eq. (A-5) Is rewritten as, 

F(top) ■ *, J >J+1 <u u+) - « 1fj ) (A-6) 

Similarly one obtains for the viscous force exerted by the stream (i,j) 
on the stream (i,j-l), call it F(bottom), 

F(bottom) * Vjj (u^j - (A-7) 

The net force on the stream (1,j) by the neighboring streams (i,j+l) and 
(1,j-l) Is thus given by, 

F(top) - F (bottom) • V^, <u 1(J+1 - u 1(j > - 

(A-8) 

By adding to Eq. (A-8) a similar expression for the net force exerted 
by the streams (i+l,j) and (i-l,j), one obtains the required relation 
for the viscous forces exerted on the stream (i,j) by the neighboring 


streams. 


APPENDIX B 


B.l Introduction 


In this appendix the listing of a computer code which calculates three- 
dimensional compressible laminar viscous flow through the lower-left quarter 
of a rectangular duct Is presented. If the flow Is not symmetric with respect 
to the center-planes of the duct, the code can be easily modified to compute 
the flow across the whole duct. The code is organized along the lines of 
the author's code [17] for two-dimensional flows. 

The code consists of MAIN program and the following six subroutines: 


Subroutine START; 

Subroutine XSAREA; 

Subroutine SIPSIM; 
Subroutine ZDZYDY; 

Subroutine VSCSTY; 
Subroutine PDSOLV; 


Assigns values to the streamline coordinates, and the 
velocity and temperature distributions, at the duct 
inlet. 

Calculates cross-sectional and surface areas of the 
streamtubes . 

Calculates mass flow rates through the streamtubes. 
Calculates the new streamline coordinates at the end of 
an integration step. 

Calculates viscosity. 

Solves a penta-diagonal set of linear algebraic 
equations. 


Subsection B.2 contains the computer listings for the MAIN program and 


the subroutines. 



B.2 Listings of the MAIN Program and the Subroutines 


The following FORTRAN names are used In the MAIN program. 


U( I ,J) * Velocity along the streamline (I,J). 

T(I,J) ■ Temperature along the streamline (I,J). 

H(I,J) ■ Enthalpy along the streamline (I,J) 

RHO(I.J) » Density along the streamline ( I »J) . 

USH(I.J) - U 2 (I,J)/2. 

Y(I); Z(J) * Coordinates of the streamline (1,0). 

DZ(J) • Z(J) - Z(J-l). 

DY(I) » Y ( I ) - Y(I-l). 

NI; NJ * Specifies the number of streamlines. Total number of streamlines 
used is equal to (NI-1) times (NJ-1), and corresponds to a 
NI x NO grid. 

Nil » NI-1 

NI2 « NI-2 

NIP * NI+1 

NJ1 = NJ-1 


NJ2 * NJ-2 


NJP * NJ+1 

NBLZ; NBLY * Specifies the number Hf streamlines assigned to the boundary 
layer at the inlet. 

DX * Si 2 e of the Integration step along the duct axis. 

SIPP(I,J), SIMP(I,J), SIMM(I,J), etc. * Mass flow rates through the various 
quadrants of the stream (I,J). The last two letters in the 
FORTRAN name SIXX are to be Interpreted as; P * +, M * 
and 0*0; thus, 

SIPM( I ,J) * 

I * J 

AX(I,J) » Cross-sectional area of the stream (I,J). 

ASJ(I), ASI(J) - Surface areas of the stream (I,J). 
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ALPHAY(J) ■ aj' ; see Section 2.1. 

Other a* d, y. end «'s, from Section 2.1, ere named similarly. 


AI(I»J) * A* j used In the finite difference equations for the velocity 
or enthalpy. 


Other coefficients appearing In the velocity and the enthalpy finite difference 
equations are named similarly. 


Other FORTRAN names appearing In the program are either defined In Section B.l 
or defined via assignment statements In the program Itself. 


Following an some explanatory notes on the various subroutines: 


Subroutine START; The subroutine assumes NI ■ NJ. The grid along J Is 

set up along the lines discussed In Ref. [17]. The 
grid along I (the lower wall) Is scaled from the grid 
along J using the duct aspect ratio, ASPR. 

Subroutines XSAREA, SIPSJM, and ZDZYDY; FORTRAN statements In these 

subroutines are almost direct translations of the 
formulae given In Appendix A. 

Subroutine VSCSTY; In Its present form, this subroutine supplies a constant 

value laminar viscosity to the MAIN program. To extend 
the code to calculate turbulent flows this subroutine 
will require a major change. 

Subroutine POSOLV; This subroutine solves a set of linear penta -diagonal 

algebraic equations by using the alternating direction 
Implicit method. The boundary conditions used are 
fa Y the function Is zero at the walls, and (b) the func- 
tion Is symmetric about the centerplanes of the duct. 


C2 


The data 1$ put in via assignment statements In the MAIN program. 

The input statements have been placed at the very beginning of the program 
and are organized into the following four groups: 

a) Dynamic and thermodynamic data; 

UCNTR * Duct centerline velocity 
TCNTR * Duct centerline temperature 
TWALLZ * Left wall temperature 
TWALLY * Lower wall temperature 
P ■ Pressure 

CP * Constant pressure specific heat 
GASK * Ratio of the specific heats 
GASR * Gas constant 
PRNDL * Prandtl number 
RENUM * Reynolds number 

DPDX * Pressure gradient guess at the duct entrance 

b) Geometric duct data; 

HEIGHT * Half duct height 
WIDTH s Half duct width 

HSLOPE , WSLOPE » Slopes of the duct walls defined as, 

WIDTH at x + AX » WIDTH at x + WSLOPE -Ax, 

HEIGHT at x+Ax * HEIGHT at x+ HSLOPE- Ax. 

XEND * Length of the duct 

c) Convergence criteria, grid size, etc. 

DPDTL = Tolerance within which the axial pressure gradient is to 
be calculated 

DXFAC * specifies the integration step; DX - (HEIGHT* WIDTH)/2*DXFAC 

MAXIT * Maximum number of iterations allowed; if this number is 
exceeded the program will stop and write out a message 

NJ, NI = Number of grid points in the cross-stream plane 

NPRNT s Governs the output frequency 




d) Data for YSCSTY and START subroutines; 

EXPU, EXPY - Exponents for the velocity distributions In the 
boundary layers along the lower wall and the left 
wall respectively. 

VZERO, TZERO ■ Optional parameters for calculating laminar viscosity. 

The program prints results at specified Intervals along the duct axis. 
The output consists of, 

• 

EX * Distance along the duct measured from the duct entrance. 

HEIGHT, WIDTH - Duct half width and half height at EX. 

ITRN ■ Number of Iterations required for the last Integration step. 

DPDX ■ Pressure gradient. 

PREC ■ Pressure recovery. 

Z(J), Y(l) ■ Coordinates of the streamline (I,J) In the cross-stream 
plane at EX. 

TAUY(I), TAUZ(J) ■ Shear stress distributions along the lower wall and 

the left wall, respectively. 

QDOTY(I), QDOTZ(J) ■ Distributions of heat flux to the walls along the 

lower wall and the left wall, respectively. 

TAUX * Average wall shear stress at x. 

QDOTX * Average wall heat flux at x. 

U(I,J) ■ Velocity along the streamline (I,J). 

T(I,J) ■ Temperature along the streamline (I,J). 

The results for U(I,J) and T(I,J) are printed in a matrix form. 



or pooa Quurr/ 


tv gi release 2.0 


MAIN 


OAfl ■ 9122* 11/34/ 29 


C 

C-—--NEMC0-U3-91*—— 

c 

C0NMCN/C0Nl/U(23,25) ,T ( 25 -ZS \ ,H( *9, 23 » ,RHOI 25.291 ,USH( 29*29) 

C CMMQN/C QM2 / 2 ( 2 5 1 , 01 1 23 1 . V 1 29 > ,0V < 29 1 
COMMON/ CQM3/ N J « NJ 1 ,NJ2 ,NJP,Nl ,N( l , N 12 ,NIP ,NBL2 ,N8LV ,OX 
CQMMQN/CON4/SIPP( 25,29! »»(MP( 29,23 ) ,SIMM( 29,25) ,5 1 PHI 25,25) .PLQftAT 
l,StPa<25.25),$lMO(25.25),SlGP<2S,25l,SlQH(23,25),StOO<25,25) 
CGMKON/CQM5/UCNTA * TCNVR, THALL2 .TWALLV , P.MNV ,V2ERO, T2ERQ* EXPV »EXPU 
CONMON/COM6/AX125,23>,ASJ123) .ASK 251 

C OMMON/CON 7/ ALPHA* 1 2 S » , BIT AY 1 25 I , OELT AY 1 29 I , GAMMA V < 23 > , AL PH A2 ( 25 1 . 
1SETAZI2S) ,aELTA2l23),CANNA2(29),VISCt(23,23> tWISCJ (25.25 » 
CCMMON/COM8/HE ISHT ,W iOTtf, M SLOPE, WSLOPE ,OL0A«ASPR 
DIMENSION All 25, 23), CM 25, 25), 88) 23, 23) ,AJ< 23,29) ,CJt 29,25 ) .00125, 
129 ) , VPJ (29 ,25) ,VP I (25 ,29 ) ,MTJ ( 29,25) ,HT I ( 29, 25) ,UU( 29,25) , 

20U< 29,29) ,UN( 25, 29), 07129,25) ,m 25,29) , 

3TAU7 (20) ,TAU2(20),CDGTVt20). 00074(20) ,HHI 29,23) 

c 

C— DATA INPUT ————————————————— 

C— DYNAMIC ANO TrtERMOOVNAMlC 0 ATA— — — — 

UCNTR-O.OOIS 

TCNTR-100.0 

T WALL 2*0.0 

TWALLY-0.0 

P-71000.0 

CP-1003.5 

CASK-1.4 

GASP >237.0 

PRN0L-G.72 

RENUM-100.0 

0 POX— 0.000073 

C GCCMETR 1C OUCT 0A7 A ——————————— 

HEIGHT-3. /8. 

WI0TH-3./4. 

ASPR-W 10TH/HE IGHT 
-SLQPS-0 .0 
HSLQPE-O. 0 
XEN0-0.5 

C--CCNVEP.GENCE CR ITER l A, GRID SI2E.ANC OTHEk (SATA——* 

OPOTL-O.Ol 
ACLMT-O. 00001 
OX FAC-5.0 
MAXIT-150 
NJ-11 
N I - 1 1 
NPRNT-l 

C— DATA FOR V I SC STY ANO SATART SUBROUTINE—————— 

V2SRO-0 .0000386 


ORIGINAL PAGE IS 

OF POOR QUALITY 


IV 01 <CIE*SE 2.0 MAIN DATE * «122* U/3A/2E (A 

TZEAQMAO.O 

EXPV*.7 

EXPU*7. 

C 

£••••• INI TI ALI ZA7I0N PCA LAMINAR (LOWS ••*•••••• 

00 20 J*l.NJ 
ICTAVI JI-2./3. 

0«L7AVU>*O.S 
GAMMAYl 01*0.9 

20 ALPNAVI Ji *2»/3. 

00 21 1*1. Nl 
ICTAtl Ii*2»/3. 

OELTA2IIIO.S 

SAMNA2(I)>0.5 

21 ALPHA2I I l*2./3. 

0EL7A2I 1 )*0.29 
BCTAZI l)*A./9. 

GAMMA2I 11*0.29 

ALPHA2U I •<►./*. 

C*** END LAMINAR PLOW INI7 IALIZA7I0N ■ •••*•*■••**•* 

C**INI f l AUIZATIONS********************************************* 1 ********* 

X AMR -2.0 
l STEP *0 
I 7RN*0 
NPl*NPRNT-l 
,\JR«NJ>1 
NJl-NJ-l 
NJ2-NJ-2 
NIP«NI*l 
NI1-NI-1 
MI2-NI-2 
OC 90 J-l.NJ 
90 C I ( 2 » J I *0 • 

OC 92 

92 C J 1 1 r2 1*0 • 

OX R> i HC 13H T »M l OTH I /2 .3 SOXR AC 
QX-OXP/912. 

OXR-OXR-ACLMT 

EX-OX 

OR •ORDX*OX 

PPLIMV-l./PRNDL 

CPINV-l./CP 

AIMVal./OASR 

SASK IN* 1./ CASK 

SNOSPO* (GASK*GASA*rCNT( >**0.9 

PK4L-P 

OP 1*0. 

OP 2*0. 


m 
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PACE E8 

OF POOR QUALITY 


i't tv si 


RELEASE 2.0 *AIN OATE > 3122* 11/34/ 23 41 

0X1-0. 

qloa-height*mioth 

CALL START 

MR IT£(6, 1*1 ) !Z(J), J-2.NJ) 

MR 1TE( 6, 1*2 1 iy(Ilft*2iN(| 

CALL XSAREAtll 
CALL VSCSTYd) 

CALL S IPS IM4Q I 
CALL 4O2V0Y 

PRF AC-0 .3 *RHC ( N I ,N J 1 *UCNT R40CNTR 
OPTOT— OP 
00 33 N-1,20 
VFJ1N.21-0. 

VFH2.N1-0. 

HTJ!N,2)-0. 

S3 MTH2.N1-0. 

C 

C*«AXl AL CALCULAT 10 NS -*•••••••*••*•••••••**••*•••• •••••••••••••*••*•••*• 

100 IST6P-IST6P*! 

EX-EX*OX 
QPTQT-OPTQT H)P 
P-P+OPTOT 
CALL VSCSTY(O) 

NPI-NPl*! 

IFCXAMP, 5T.il NPI-NPRNT 

1FINPI.lt. NPRNT) CO TO 200 

NPt»0 

OPOX-QP/DX 

PREC-OPTQT/PRFAC 

T AU-0 » 

QDOT-O. 

00 l Li 1*2 ,Nt 
TAUVm-2.*VFJ( l»2)«UU,2) 

UD0TYm-2.*HTJ< !,2>*>(H! 1,21-HU, II 1 
TAU-TAU*TAUY! (> 
m aoor-gooTtwooTYm 
OC 112 J-2 ,NJ 
TAUZUJ-2.*VPl !2»J)*U<2,J) 
iQGTZ! J )*2 .*HT l(2,JJ*(H(2»JI“H! L. J » 1 
TAU-TAU+TAUil J) 

112 ODOT-QDOT'-OCQTZU) 

TAuX-TAU/( HEICNT+W1CTH1/0X 

S0GTX-fl03T/<H6 IGHT«* 10 TH I /OX 

*iFlTE(6, 1301 EX »H6 10HT , wiO TH, 1 TRN *OPOX ,PRF,C 

130 FCRMAT(lHO,'OIS. ALONG X-',F7.*,'*» HEIGHT-' , FT. *,'••• WIDTH-', 

1 F7. *.'*** ITERATIONS-' ,13, '••■OP OX-' ,611 .*, • — PREC- » , 6U.*1 

..aired, i3i) taux.qoqtx 

131 FORMAT! 1H ♦ , TAUX»',611.*,'*»*QOOTX«',6ll.*> 



I 



>< I'/ SI RELEASE 2.3 "ft!* 


OAVS • 91224 U/36/29 »i 


*AITE(*.14l) (ZUI.J-2.NJ) 

•A !TE( 9.1421 (YUI.l-2.NU 
•rite(*.i44) (tauyu i,i-2,nii 
-MTE<*.149) (TAUZUI .J-2.NJ) 

PRITE(*»14*> ! GOOTY! I). I -2 ,Nl) 

S.*lTSi*.l47l (OOOTZ! Jl.J-2.NJ) 

DC IDS JREV-l.NJl 
J-NJP-jREV 

139 ..AIT6(6.14JI J,(U U.JI.I-2.NU 
DO 13* JREV-i.NJt 
J-NJP-JHCV 

136 WAI TE!*.149) J.I T< I.J ) .1-2. NI ) 

141 FORMAT ( 1HO « 1 Z* ,)|X,lPtOEll.4/l 10X, IP 106 11.4) I 

142 FORMAT ( 1H , • Y* . 0X. IPlOEll .4/ < 1 OX. 1 P 10E 11 .4 ) ) 

143 FORMAT ( IH ,'li J-* .12.2X. tPlOEli.4/t 10X . IP 10611 .4) ) 

144 FORM ATI IH . * TAUV ,3X, lPlOf 11 . 4/ ( 10X , IP 10611 .4 ) ) 

. 149 FORMATUH . 'T»UZ'.9X,lP106ll.4/( 10X.IP10611.4I) 

146 FORMATUH ,»OOOTY*,4X.lPlOEll.4/( 10X , IPlOEl l .4) ) 

147 FORMATUH , ‘80072 • »4X, IPIOE 1 1.4/1 10X , IP 1061 1 .4) > 

140 FORMAT! H ,'T J-‘ . 12 . 2X, 1P10E1 1 .4/ ! IQX , IP 101 U .4) > 

(F (EX.GT.XENO) GO TO 3000 

C 

c— resetting dynamic. thermo, ano GEOM. variables for next iteration 

200 ITRN-t 
X 1-0X1 
X 2- OX 40X1 
SUMX-X14X2 
SUMXX-X l-X l4X2*X2 
SUMP'-OP 4OPt*0P2 
SUMXP-Xl-OP 14X2*0P 
0EN-3.«SUMXX-SUMX«SUMX 
JPSL0P-! 3 .-SUMXP-SUMP- SUMX I /QEN 
DPZ6RC-(SUMP-SUMXX"SUMX*SUMXP ) /OEN 
0X1 -OX 
OP 2-OP l 
OPl-OP 

IF! 1STEP.LE.2) GO TO 220 
IF(QX.GE.OXF) XAMP-l.O 
UX-0X-XA4P 

1FUST6P.LT. 4) GO TO 220 
OP-OPSLUP-IX2 40X) *0PZ6R0 
220 PSTART-P 

OPOAMP-O. 9 -EXP (-0. 115-6X/H6 IGHT) 

UOLO-U! 2. 2 ) 

CALL SIPSIMIO) 

CALL XSAR6A10) 

ZY FAC -FLORA T-U( NI ,NJ ) / (CIOA-FLORAT-W N t ,NJ I -GASK l.V/P ) 

ZY TLRN-OP l-OPOTL-OLO A/ ZYF AC 


m 
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ORIGHv'U 




OF POOR QUALITY, 


*.N .y 61 ft Cl CAS 6 2.0 NA IN OATS • 8122* tl/M/21 

UTUN*OPl*OPOTl*etDA/.PLQRAT 

IP (2VTIAN.IT »ACIMT*OLOA) ZYTlAN-»e!„MT*0l0A 

ip (uruN.iT,Acifir*uoioi utirn-aciht*ucid 
c 

C**** EVACUATION OP R.H.S. OP MON* ANO ENERGY EOS. •••••*•••»••••••••••• 

00 300 4-2, N4 
00 300, 1-ZrNI 

VP4U,4)-0.J*AS4(I)*VISC4I l,4)/0Z(J) 

VPII I,4)-0.5*ASK4)*VI$C1( l,J)/OY( I) 

HTI(l,4)-VPt< l,4)*PRUNV 

300 HT4(t,4l*VP4( 1,4)*P«UNV 
OC 301 4-2 ,N4 

NT l (MtP,4)-HT KNI.4) 

301 VFltNIP,4)-VFI(Nl,4) 

00 302 l-2.NI 
HT4(I,N4P)*HT4U,N4) 

302 VP4Il,N4P)-VF4< l,N4) 

0U(2,2 )-U( 3.2 ) »(SIPP(Z»2) ,*AIPHAZ( 2)*SIPM(2»2)*VF 1 1 3* 2) )♦ UI2.21* 
U2.*SlPPI2.2)*3.*AtPHAV(2)*SINPI2,2)*3.*AlPMAZ(2)*SlPH(2,2)* 
2AlPHAZm*SIMH<2,2)-VP4(2.2)-VPl(2,2)-VP4I2,3l-VFIl3,2l) 

3 ♦Ul2i3)P(SIPP( 2, 2)*AIPHAY(2)*SIHP( 2.2 ) ♦VPUI 2.31 ) 

OO 310 1-3, NI 

310 OUI I ,2) -U( I~l ,2 ) *( SIMP! 1 , 2) *AIPH AZ ( I) *5 INNI 1 .21 ♦VP It 1,2) )♦ 

1 UIl*l»2)*(SIPP( 1 ,21*AIPHAZ( I 1 *SIPN< I. 2) ♦VP I ( !*l,2) I* 

2 UIl, 2) *12. *51 OP (I, 2) *3.*ALPHAZ( I J *SIOM( 1 ,2 )«*VPI ( I ♦1.21 "VP J (1*31- 

3 VPK I.2I-VFJI 1.21 IHII I, 31*151 QPt 1.2) ♦VP J( 1.3)1 
OC 311 4-1, N4 

311 0U(2,4)-U(3»4)*(SlPO(2»4)*VPt(3»4))H)(2»4)*(2.*SlPQ(2»41* 

1 3.*AIPHAY(4)*SIM0(2,4)-VPI (2i Jl-VPJI 2,4 ♦) )“VPt (2.41-VP4 1 2,4 1 )♦ 

2 U(2,4*1)*(S IPPI 2 .4 1 ♦ALPHA Yd *, ^51 MP ( 2 .41 ♦VP41 2.4*1 ) )♦ 

3 U<2,4-l)*(AlPHAV(4)*SlMN(2,4)*StPM(2,4)*VF4t2.41) 

00 312 4*3. N4 

00 312 1-3. NI 

312 DU( 1 .4 )-U( I* 1 . 4)* ( SI POI 1 .4 1 ♦VP I IX>1 .4) I ♦UI I- 1 .4 )•( S IMOI 1 ,4 )♦ 

1 VPK 1.4) (♦UI I,4)*(SI00(l»4)-VFI( 1*1,4 )-VF4( 1,4*1 )-VPI ( l .4)- 

2 VP4II.4) )HI( t * 4*1 )*( SIOPII »4)*VP4( 1 ,4*1 ) 1 HI 1 1 ,4-1 1 *( SIOMI 1.4 )♦ 

3 VP4U.4I1 

00 315 4-2. N4 
DC 315 1-2, Nt 

315 00 ( 1 .41 -OU ( 1 , 4 )— QP*AX (1.4) 

C 

QT(2»2)-H(3»2)*(SIPP(2»2)*B6TAZ(2)*SIPN(2»2) *HT I (3,2) ) *H(2,2)* 
1(2.*SIPP( 2,2) *3.*B6TAY(2)*SIMP(2,2)*aETAZI l)*SIMN<2,2 )*3 . *8ETAZ { 2 ) 
2*5 IPN( 2, 2 ) -HT I ( 3.2 )-*HTJt2» 3 )-HT l (2, 21-HTJ I 2,21) ♦HI 2, 3 ) *•< SIPPt 2 .2) 
!♦ SET AY( 2) *S IMP( 2 , 2 )♦« TJ( 2 , 3) ) *USH( 3,2J*(SIPP(2,2)*0EITAZ<2)- 
*SIr"1(2,2) ♦VFK 3.2) ) ♦USHI 2 , 2) *(2 . *S IPP ( 2 . 21 ♦3.*0ELTAY (2 ) *SIMP(2 , 21 
5 ♦OElTAZt L)«SI.MH< 2,2) ♦3.«0£LTAZ(2 ) -SIPM( 2 ,2 )-VFH 3 ,C )-VF J I 2, 3 1 
a -VP I(2,2)-VFJ(2,2)I ♦USH(2,3)*( S IPP( 2, 2 i *06tTAY(2 ) *SINP( 2, 2 ) 


ORIGINAL 'V:*; £0 

OF POOR 


AN IV 01 ACLEASE 2.0 MAIN DATE *91224 U/3A/2P 

7 *VF0(2.J»> *H( l»2)*HTt(2.2)*H(2«ll*HTJ(2.2) 

OC 321 I>3.NI 

321 3Tll.2)*M(l*1.2)*ISIPP(I.2>*8«TA2ll)*SIPM(l.2>*HTt(l*l,2n * 

* lHU«1.2)*IStMP( l.2)*ftETA2t II*$tMM(I,2)*HTI< 1,2) M> H(t,2)*<2.* 

2SI0PI 1. 2)*3.*8STA2(1)*SIQN<1 ,21-HTl (1*1. 2)-HTJ(l,3 l-HTI < 1.2)- 
3MTJII,2)I *H< I,1I*ISI0P(1,2)*HTJ<1,3)) *USHI 1*1 .2 )*(S IPPI 1.2)* 
40CLTAZI !>«SIPMI t.2)*VPIt l*l.2>> *USHl 1-1 .2) • ( SIMP! 1 .2 >*OELTAZ( I ) • 
SSIMMI 1.2) *VPI( 1 . 2 ) ) *USH(I #2)*(2.*SI0P( 1.2I*3.*0EITAZ(I )*SION( 1 .2) 
A-VFl(I*l,2)-VFJ< I.3I-VPII 1 .2J-VPJ ( 1 .21 1- *USH|l,3)*(SlOP(I.2)* 
7VP3II.3I) «H(t.U«HU(!.2> 

OC 322 0*3 .NO 

322 OT(2.0)*HI 3,0 > *1 SI PCI 2 .0 ) *HTl (3.0 II *H( 2.01*12. •SIP0I2.J 1*3. • 
18CTAVIO)*SlMOi2 »0)-HT 1(3.0 l-HTOt 2.0 *1)-HT! I 2 .OI-KTOI 2 .0) ) * 

2HI 2.0*1 )• (StPP(2»0>*SETAV< 0)*SlNP(2,OI*HTU12.0*ll I* HI2.0-1)* 

3 1 SET AVI 0)*S INHI 2* JI*SIPN(2 .0) *HTO 1 2.01 I* USHI 3 , J> *i SI POI2 .J ) * 
4VPII3.01I* OSH 1 2.01*1 2 .*SIPQI2 .0 ) *3.*0ELTAV( 0 )*SIMQ(2.0 )- 
9 VP 1 1 3 .0)- VPOI 2 .0*1)-VP( 1 2.0I-VF0 1 2.0 1 ) * USHI 2 .0*1 > *1 S IPPI 2 .0 1 * 
60ELTAVI 0l*SIMP(2»UI* VP 0(2. 0*111* USHI 2 .0-1 > •IDELTAVIO ) • 

7S I MM 1 2.0 ) *S IPNI 2.0 l*VPOI 2.0))*' HI 1.0) *HT I (2.0) 

00 323 0*3. NO 
00 323 1*3. NI 

323 OTI1. 01 >HI l*l,0>*(SIPC(t.O>*H7I(I*l,01)* HI l-l.O)*(SlNail.O)* 

IHTI 1 1 .0))* HI I .01*1 SI OQI l »0)-HTI I t*l.OI-HTO< I .0*1 )-HTI 1 1 .0 )- 
2HT0I 1. 0))* HI 1 .0*11*1$ IOPI 1.01 ♦HTJI 1.0*1))* HI 1 .0-1) *( SIOHI 1 .0) * 
3870(1.0))* USHI I*1.0)*(SIPOII ,0) *VPII 1*1 .0 )) * U$H(I-l.OI*( 

4SIMCI 1.0) *VPI I 1.0)1* USHII,0)*(SIOOII.O)-VPIII*l.O)-VPJ(I,0*l)- 
5VPIU .OI-VPOI 1.0) )♦ USHI I .0*1 )*( S IOPI 1 .0 )*VPO( 1 .0*1 ) ) ♦ USHI 1. 0-11* 
6ISI0MI I.OI*VPO( 1.0)1 

C**l TEA ATI VC CALCULATI CNS POP NEW U .H .2 ••••••••••••»•••*•*•••••••*••* 

C**CALCULATIQNS OP NEo VELOC IT US*************************** 

00 70 900 
400 CON 7 INUE 

00 410 0*2. NO 
00 410 1*2. NI 

VPOI I .0) *0 .9*4 SO I ( )«V I SCO (I.0)/C2IJ) 
VPltI.O)*0.9*ASI(0)*VISCl<I.O)/OVm 
HTK I.O)*VPI( 1 .0I*PPIINV 

410 HTJl leO)-VFJI l.OI-PRLINV 
00 411 0*2 .NJ 
HTHNIP.J) -HTIINI .0) 

411 VP 1 1 N IP »0)*VPH NI .0) 

00 412 1*2.. NI 

HfJI I,NOP)*HTJI I. NO) 

412 VPOI I. NOP) -VPOI I, NO) 

900 CONTINUE 

AU2.2I— SIPP(2.2)-ALPMAZ(2 )*SIPM 1 2 ,2 ) *VPI I 3 .2 ) 

38(2.2) -2 .-SIPPI 2* 2) *3 •* ALPHA VI 2) *S IMP 1 2.2 ) *3. -ALPHA! ( 2 I -S I PHI 2 .2) 
l* ALPHA! I l)*Sli y *(2»2) *VP J (2.2)*VFI(2.2) »VPJ 1 2 ,3 ) *VFI (3.2) 



60 / 

ORIGINAL PAGE !S 
OF POOR QUALITY 


(V 41 RELEASE 2.0 MAIN 9ATE • 9122* 11/16/29 Pi 

A4(2.2)* -S IPP( 2 * 2 l-ALPH AY ( 2) *SIHP( 2,2 1 *VP J< 2, J 1 

00 915 I - 3 , M l 

*1 (I ,21 ■ -SIPP( t,2)-ALPHAZ(I)*$lPH( t ,2 )*VP III * t ,2 ) 

CIII.2I* -SIHPI l,2l-ALPHAZ(ll*SIHH< 1,21 *VF I ( 1,21 
99(1 ,2I*2.*$I0P( l,2>*3.*ALPHA2Ci)*StON< U2 I *VF I < !♦! ,2 )♦ 
l '/P4( t,3)4VFI( t,2)*VFJ( 1,21 
A4 ( 1 , 2) ■- S 10P( 1 ,2 )*VP4 11,3) 

515 CONTINUE 

DC 516 4*3, N4. 

Al(2,4)*-SIPa(2,4>4VPI (3,41 

a9!2,4)*2.*SIPCI2,4)*3.*ALPHAV(J)*$ IMO ( 2,4 1 *VF t ( 3,4 ) ♦ 

1 VFJI2,44l>4VFU2.4»*VFJ<2,4> 

AJ (2, 4 )■ -SIPPI2,JI-ALPHAY(4)*SlNP(2,J)*VF4I2,4*l> 

CJ < 2 , J } « -ALPHA Y< J)*SIMM(2»4)“SIPH( 2, J (♦VP J( 2, J ) 

515 CONTINUE 

00 520 J- 3,NJ 
00 520 1-3 , Ml 

A( It ,4)*-SIP0< t,4M>VFl (1*1,4) 

CII l,4)*-SI*0( l, J I ♦VP l ( 1,4) 

99 1 l,4)*S I00( I ,4)*VPt( l*l,4)*VPJU#4*l)*VFI( t,4)*VFJ(I,4) 

A4( I ,4)*-S ICP( l»J)*VFJ<I,4*l) 

C4(l ,4»*-StOM(t,4l*VP4< 1.4) 

520 CONTINUE 

CALL POSQLVI Al,Cl,8B,A4*CJ ,00,U,UU) 

IF ( ITRN.GT.l) CO TO 340 
00 530 J-2.NJ 
00 530 1*2, Nl 
530 U(1 ,41 -UU( 1,4) 

540 JO 550 4*2, NJ 
00 550 1*2, NI 

550 U(l,4)*0.5*(U( I * J) *UU( 1,4)1 
00 560 J*2 ,N4 
560 U( N iP , 4 ) "U (NI 1 ,4 ) 

DO 570 t*2,NI 
570 J( l , Y4P) *U ( I ,N4l ) 

U (Nl 0 ,NJP ) *U( N I l.NJl) 

C 

C**CALCUL ATIGN OF NEW ENTHALPIES**************************** 

00 610 J-2.NJP 
DO olO 1*2, NIP 

olO USH( t,4)*0.5*U( I , 4 )*Ul t,4 ) 

A[(2,2)«- SIPP( 2, 2 ) -9ETAZ ( 2 ) *Sl PN( 2 » 2 ) *HTI (3,2) 

8B(2,2)*2.*S1PP(2«2)*3. *96 T AV ( 2) *S IMP ( 2, 2 ) ♦SETA 2 ( l ) *S INM( 2 ,2) ♦ 

13 . *9£T A2 ( 2 ) *S 1 P“( 2,2 )*HTI (3,2) ♦HTJ ( 2, 3) *HT I ( 2 , 2 ) ♦HTJt 2 ,2 ) 

A4 (2 ,2) * -S IPP ( 2 ,2 )-9ETAY( 2 ) *S IMP ( 2 , 2) ♦HTJI 2,3) 

TT( 2,2)*0T(2,2I *H( i,2)*HTt U^l^HUaJ-HTUU^HlSHO^m- 
lStPP) 2,2)-06LTA2(2)*SlPM(2,2)*VPI( 3,2) I -USH( 2, 2) *(2. *SIPP(2, 2) ♦ 
22.*QELrAY{2)»SIHP(2,2)H36LTA2(l)*SIMN(2,2l^3.*OELTA2(2)*StP“(2,2)* 



6 


ORIGINAL M/Uje ;<i 

OF POOR QUALITY 




(V 31 ML64S6 2.3 MAIN DATE ■ 5122* 11/36/24 »l 

3VFH3.2)*VFJ<2.3)*VPI(2,2)*VFJ<2.2) )* U3H( 2 , 3 ) -I-S IPPI 2 .2 >- 
•*06tTAV(2>»SI^(2,2>*VFJ<2.2n 
00 619 t-3.NI 

At (1, 21— SI PM l,2)-9ETAZ< | l-SIPNI 1 .2* -HT 1 1 1*1,2) 

Cl It ,2)— SINPI I,2)-tfTA2U I-SINNI 1,2,) *HTl( 1.2) 

9B(1 ,2 )-2 . -S I Cat t »2)*3.*9ET A2( D-SIONI 1.2) *H Til l*l,2)*HTJ(l,J)* 

lHTI(I, 2 )*Hrj( (. 2 ) 

AJll ,21— SIOPU,2l*NTJU,3» 

TTt I,2l-OTU,2l*HII»l)-HTJ(I,2)* USH(l*l,2)*(~S(PP(I , 2I-0ELTAZ( ()• 

ISIPK l,2)*VFl< 1*1,2) ) * USNU-t,2)-(-S INPU .2)-06LTA2U) -91**11 ,2) ♦ 

2VPIU.2II- JSH( t,2)-(2.-SI0P(l ,2J*3.*0IITA2( 1 ) -S I OP,< 1,2)* 

3VF l ( 1*1,2 )*VPJ ( 1 , 3 )*VPU I , 2)*VPJ( 1,2) )*USH< 1 ,3)-(-SIOP( t ,2)* 

4VP J( 1,3)1 

619 CONTINUE 

00 616 J-3.NJ 

AI (2,JI— S 1P0I 2.J)*HT1(3.J) 

eSI2,J)>2.-$ tPOI2,J)*3.-6CTAV(J)-SIMO(2,J)*Hrr<3.J)*Hrj(2,J*ll* 

ihtk2,j)*htj(2,j) 

A412.J)— SIPP(2,J)H»fTAYU)-SlNP< 2, J) *HTJ( 2, J*1 ) 

CJ (2 »J) — 9ETAV ( J )-SINM(2, J )~SI P*( 2 #J) *HTJ( 2. J ) 

TT(2,JI-0T(2.J)*H< 1,J)*HTII2,J)* USM(3.J)-(-SIP0( 2, J ) *VF 1 1 3. J ) )- 
1USH( 2,J)*(2.*SIPO(2,J )*3. -OCLTAV (J)-SIM0(2,J)*VPl(3,J)*VPJ(2,J*l)* 
2VPt(2,J)*VFJ(2,J) )* USH(2,J*l)-(-SIPP(2,JI-0fLTAYCJ)-SINP«2,J)* 

JVFJI2,J*l> )♦ USHl 2 ,(J“l )•( -OELT AY I J ) - S IMM( 2,JHSIPM(2,J)* 

4VFJC2,4)) 

616 CONTINUE 

00 620 J-2 ,N4 
00 620 1-2, NI 

AI(I,J) — SIPOI t»J)*HTI (1*1 ,JI 
CHl.Jl — SlNOII.J)*HTl <I,JI 

99(1 , J )-SIQQ( I , J) *HT I ( I*l,J)*HTJ(l,J*l)*MTK l,J)*HTJI 1,3) 

AJ( I*J) — S IOP( t »JI*HTJ( 1,3*1) 

* Cjtl.J) — SlONU,J)*HTJtl,J) 

TT(I ,J)-OT( I , J ) ♦ USH< I*l.J»-(-StPO( r,J)*VFt( i*i , j) > * ush ii -i , jJ- 
n-siMuii , j)*vfi ( i,j) )-ush( i ,j)*< sion< t ,j)*vfk i*i,ji*vfji t ,j*u* 

2VFM I ,J)*VFJI I.J) l*USHII,J*l)-(-S IOP( I,J)*VFJ( t , J*l ) ) *USH U . J- 1) • 

3(-ilCN( I,4)*VFJI t.J) ) 

620 CONTINUE 

00 621 J-2.NJP 
00 t>21 1-2. NIP 

621 HU(J)»HII ,J)-Mi,i) 

CAui. posolvui , ct ,ao,Aj,cj,rr.H.HH) 

IF < ITPN.OT.l ) 30 TO 640 
00 630 J-2 »N J 
00 630 I -2 ,Nt 
630 M ( I I.J ) 

040 00 o 90 J-2.NJ 
30 690 1-2, Nt 



62 


OTOGIMM- 
OF POOR QUA-t ‘ * 




AN tV 31 RELEASE 2.(3 MAIN 3ATC > *122* 11/36/2* »< 


650 ( t,J»-0.5*(H« l ,JI*HH< I.J) >*H< 1,1) 

00 6*0 J-2.NJ 
6*0 H(NIP,J)-H«NI1,J) 

00 670 1-2, Ni 
670 HU.NJU-Hl I.NJl) 

HI NIP ,N JP ) *HI M l ,NJI ) 

00 U00 J-2.NJP 
00 UOO 1-2, NIP 

uoo rn,ji*M( t, j)*cpinv 
c 

C«* CONVERGENCE CHECRS*****************************" 1 ************** - ***** 

CALL ZDZYDY 
ZCQR a Z (NO) -HEIGHT 
YCOR-YINIHWIOTH 
UCHEK*U( 2,2 l-UCLO 
UCHEK-AISt JCHEK) 

Z YCHEK«A*$ I HI OTH*ZCOR ) *A*S (HEIGHT »YCOR ) 

I PIUCHEK.GT .UTLRN1G0 TO 950 

IP ( ZYCHIK. LT .ZYTLRNI GO TO 100 

UPCClfU-<ZCa**W!OTHtVCOMHtIGH?)*ZYP*C/OLOA 

OPCOR *OPC OR *0 P0AHP 

QP-QPACPCOR 

00 9*0 J«2,NJ 

00 9*0 l-2.NI 

9*0 00lt,JI«0Ul I,J)-0P9*XII.JI 
950 IPIITRN.ST.HARIT) GO TO 2100 
U0L0-LI(2,2) 

CALL VSCSTYIO) 

ITRN*lTRN*l 
30 TO *00 
2100 UR 1TE(*,2101 ) 

’ 2101 PORHATI/ » IX ,' ••VELOCITY-PRESSURE ITERATIONS EXCEEDED THE LIMIT SET 

1 8Y MAX IT | 

3000 CCNTINUE 

STOP 

EHO 
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OF PO&.S Q1./1LV 


RELEASE 2.0 START OATf > S1224 11/34/29 Pi 

subroutine start 

CSMMCN/C0Ml/U(2S,29l,T(23,29l »H( 29.29) .RHC129.29I »USH<23,29» 

COMMON/ CQM2/ Z ( 29 ) .1)2(29 1 . V (29 ) .CV ( 29 ) 
COMNON/CaMS/NO.NOl.NOZ.NOP.NI.NIl.NU.NIP.NRLZ.NILV.QX 
C0MM0N/CQN9/UCNTR.TCNTR.TWALLZ.TWXLLY.P. AINV.VZERO.TZERC.EXPV.EXPU 
C3MMON/C0M6/ A X( 29 .29 ) . ASJ ( 29 1 . AS I ( 29) 

COhMCN/CGMA/ HEIGHT .4 I0TH.HSL3PE.mSL CPE .0L9A. ASPR 

CP-1003. 9 

ISL-O.l 

VI SKIN*VZERO*l 300.0/TZER0»M»EXPV/ (P*R INV/300 .0 I 

PNO*PLGAT«NOJ 

U2Z2*L200.*V( SKIN 

00 LO 1*1. NIP 

RHQ( I . l)*l.O 

Tii,u<rmi 

H ( I » l J *TMA LLZ*C P 
USH( I . 1) *0.0 
10 UII.U*0.0 
00 19 0*1. NOP 
PH0(l,0)*l.0 
T( l,GI«TMALLY 
M( l.OI»TWALLY*CP 
USH( l. 0)«0.0 
19 U(1.JI«0.0 
ll l)*0.0 

• ozm»o.o 

UZ2*»UCNTR**6XPU*U2Z2/Z9L)M»( l ./U.*EXPUl ) 

U22*( 10.00**€XRU*U2Z2/ZBU *•(!./« l.*EXPUH 
Z(2>* U2Z2/UZ2 
0Z(2)*Z«2)-ZI 1) 

EXPO* ALOO( HE IGHT/ Z (2 1-1.1/ ALOGC PNO/ 3. 1 

00 20 0*3. NO 

P0*PLGAT(0» 

£( J)*Z(2>*(HEIGHT-Z(2) )*(PO/PNO)**SXPG 
23 OZ(0)*Z(OI-Z(O-lt 
UC 30 I *1 #NI 
r( n*zm*ASPR 
30 OYUI*OZI I l*ASPR 
00 ISO 1*2, NIP 
00 190 0*2 .NOP 
BHOU. 01*1.0 
TU.OI-TCNTR 
H ( I ,OI*TCNTR*C? 

USH( 1 .0 1*0.9 *CICNTR*OCNTR 
190 0(1.01 *UCNTR 
RETURN 
ENO 
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AN IV Cl A CL CASE 2.0 XSARfA DATE • *1224 U/J4/2* 

SU9R0UT INK XSAACAINR) 

COMMON/ CON l/UI 29.29) . T( 29 .29) .HI 29 .29) ,AHO< 29.291 ,JSH(29, 29) 
CQMMCN/COM2/2I 29) .02(29). Y( 29) »0V 129) 

C0MMCN/CCHJ/NJ,NJI,NJ2.N.)F.NI,NII ,Ni2,NIP.NIL2.M*LY.0X 
COMMON/ C0M4/ AX ( 29.29 ) . ASJ (25) .AS I <291 
CCMNON/COMi/Nf ICHT 'U (0 TH. N SLORC 'MSLflRC'OLOA ( ASM 
C 

C« CALCULATION UR NIW ‘OUST CROSS-SCCTION •••••••••••••••••• 

HilOWr-MilGMr+HSLORf^OX 

MlOTH-NlDTH+WSLORC+OX 

olca-mcichtawidth 

C-- SUMP ACC AREAS FOR OX-l 

4SJ(2)-(V< JI+VI2)>*0.9 
00 10 I >3 'Nil 

10 ASJUMYU+U-YI l-l )) *0.9 
ASJ ( NI )•( Y (NI )-Y( Nil ) ) 

AS I ( 2)-< 2 ( J) +2(2 ) >"0.9 
00 20 J-3 iNJl 

20 AS! (J>-(2( J+ll-21 J*l) )*0.9 
AS l(NJ)-< 2 ( NJ )-2(NJl i ) 

C"» CROSS-SECT IONAL AREAS *•»••••••••••••••••••••• 

OC 90 J-2.NJ 
00 90 1-2. NI 

90 AX ( I ' J) -AS J( I ) "A S I ( J ) 

C*» SURFACE AREAS •••••••*••••••••••••••••*••••••* 

00 *0 I-2.Nl 
60 a,;j( D- 4 S jm-ox 
OC 69 J-2.NJ 
69 AS 1 1 J)-AS I ( J ) "OX 
ME TURN 
CNO 


ORIGINAL PAGE IS 
OF POOR QUALITY 


VI IV „l RELEASE 2.0 SIBSIH OATI • BUZA 11/3B/2I 

SJBBCUTINE SIBSINI'IBI 

CGNNCN/C3NI/UI 29*29) »T( 29, 29), HI 29, 29) ,BMQ| 29*291 ,USH<29,29) 
CCNNON/C3M2/ZI 29 1 *021291 *V 1291 *0V ( 291 
CUNHON/CQNJ/NJ,N.li,NJa.NJN,NI,Nll,NI2,N|N,NOLZ,NOLY,OX 
CQNMON/CON4/SIPN<2S,29>»SiHB(29,29l,SlMN<29,29>,SINN<29,29),BLQNAT 
USING! 29.29). SINCI 29*291* SJONI 29. 29) *S ICNI 29,29) .SICOI29, 291 
CQNNCN/CON9/UCNTN,TCNTN,T-ALL2,TWALLY,N,NlNV,V2ENO,T2INQ,2i<NV,EXNU 
CONNON/CON*/ AX 1 29. 29 ) . ASJ ( 29) , AS I ( 29) 

CONNON/COH7/ ALPHA Y|29 ) * IKT AV 1 29 ) , OELTAV ( 29 1 * CANNA Y (29 ) , ACNNAZI 29) * 
IICTA2I29) .OELtAZI 29), GAMNA2<29),Vt SC! (29*291 ,V ISCJI29.29 ) 
CQNNCN/CCNB/MKGWT *J10TH»HSLQNE .WSLONE »OL0A* ASNB 
3 (MINS ION BUI 29*29) 

D2INJNI *021 NJ) 

UYININI-OYINI) 

NNB-l./AA. 

H022-02I 2 )/32 • 

H0Y2-0VI 21/32* 

C-- SI'S* EXCEPT SINMI 2,2) , NEOUCED BY A TO ABSOBB 1/A ABF CAB INC IN NOE- 
30 .10 J-2.NJN 
00 10 02, NIB 
10 BUII *J)-BHO( I*JI*U(I*J) 

S!NN<2*2)-GANNAZ<LI-NUI2,2)*Y(2)-2(2I 
SINN(2,A)-<2.*NUI2,2)*NU(J.2)*NUI2,3)»*OY|3l*OZ(3l*NNN 
SiMNI2»2)-UANNAY»2l*(3.*NUI2.2l*NU(2*3)l*H0V2*0Z(3l 
SINN(2,2I-6ANNAZ( 2)* I 3 ••BUI 2 ,2 I tRUI 3,2 ) )*0V( 31 -HO 22 
BLCRAT-S(NH(2,2 )»SI BN 1 2,2 )+SINN 1 2,2) *S INN( 2,2 1*0.29 
00 20 I>3,NI 

SIBBI 1 ,2)-(2 « *BU( 1,2) *BU( 1*1* 2 l*NU( 1,3) ) *0V| I *1 )*02( 31 *BBB 
SIMB(I,2)»I2.*BU( 1,21 *BU( 1,3) ♦RUI 1-1,2) )*0V( I ) *0213 )*BBB 
SI BN 1 1 ,2) -GAHNA2I I )•( 3 .*BU U,2)*BU( 1-1,2) l-OVI I ) *H0Z2 
SI BN 1 1, 2) -0 ANN A2 1 I )• (3 .*RU ( I .2 >*BU( l *1 ,2 ) ) *0V( l*l)*M022 
20 NLCRAT-NL0RAT*SINNII,2)*SIHN( I, 21 tS INNI I,2)*SINN( I ,2) 
NLORAT-NLORAT-SINNINI ,2) -SI BN (M ,2) 

DU 30 J-3.NJ 

SIBB(2,JI*(2,*BU(2,J)*BU(3*JI*BUI2.Jtl))*OY(3)«OZU*L)*BBB 
S IMB I 2, JI-GANNAV I J |*(3 **RU( 2, J ) *RU( 2, J*l ) ) -HOY2-OZI4* l ) 

SINN ( 2, J) •OANNAYI J )•( 3.-RU I2»J)*RUI2»J-1) )*H0V2*02 I J ) 
SIBNI2,J)"(2.*BU(2,J)*BU(2,J-U*BU(3,J))«0Y(3)*02( J)-NBB 
30 BLOB AT -FLORA 7 -SI BN (2, JI-SINN 12, Jl -S INNI 2 ,J I -S I BNI2, J ) 

BLOB AT -BLOB AT-S INN( 2,NJ )-SIBB( 2 ,N J) 

00 AO J-3 ,NJ 
00 39 I-3.NI 

SIBBI I ,JI-I2.-RUI l,J)*BUII*l,JI*BU( I*Jtl))«OZ( J-l ) -0 V< I -l I -BBB 
SINFII.JI-I2.-BUI I.JI-ftUI I-1,J)*RU<I , J- l) ) -021 J-l >*0Y( l ) -BFB 
SINN I I , J ) ■ (2 • -BUI I , J ) *RU« I— l . J ) *BU 1 1 , J— 1 ) ) *021 J )-0Y|| I-BBB 
SIBNi I,J)-I2.-BUII,J)*«U< l*l,J)*BUII,J-l))*02IJ ) «0 Y ( I ♦ l ) -BBB 

39 BtOKAT -BLOB AT-S IBBI l.J USINNI I »J ) -S INNI l»J ) -S I BNI I , J ) 

40 PLCXAT-BL0RAT-SIB*INI ,J)-SIBNINI.J) 
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OF PCuU Vf 


MLEASE i.C SI PS 1.1 QATK • 91224 

CU 90 (OtNI 

90 Pl0HAr»PU3BAT-StnPt t ,NJ)-9tPPt t *.M4» 

PIOAAT-PIOAAT ♦StPPINt ,NJ> 

PU}AAT«Pl,0RAT*4. 

00 99 4«2>NJ 
OC 99 l *2. Nl 

SIP0U.4)-SIPP(I,J»*SIP*I1,J» 

SIIQI I,J)-S11PU,JI*SIM"<I ,41 
StCP(l,4l-StPPt t f 4»*S!NPtl ,41 
StOMIt,4)-SI,1*< I,4)*SIPHU,4) 

99 SIC,QII'4>-2.«<SIP0<t'4l»SIMailt4>) 

tptip.n.i) sc ro too 

WA!T6IA,43» PCCAAT 
A3 POANATUM , ••••PU3AAT-* , IPtll .41 
00 70 4AIV-l,N41 
4»N4P-4MV 

WAirE(6,69)J,(StHP((fJ),SIPP<I,J),I«2iNI) 

70 .<UT!(6,*5IJ,ISI-1(l,J),5!PH< t.J) , I-2.Nl) 

A3 PGAMATt IH J»Sl2,3X,lP10«ll.4/< UX , 1P10C11 .41 > 
ldO AfTUAN 
(NO 


11 / 34/21 »< 


t 
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V 01 KSUiASE 2.0 2D2YPV OATS ■ 91224 11/36/ 21 >< 

SUtttOUTIME 202 YOY 

CCMQN/C0Ml/U(29,23>»T(23.29),H(2S.23),*NQ(29.29»,USHI23,2S) 

C0NMCN/C0H2/ 2 ( 29 1 , 02 1 29 ) , Y ( 23 > ,0Y < 23 ) 

CCMMON/CQMl/NJ.AJl.NJa.NJP.NI .NI l ,N (2 »?) IP.NOU.NOCY.O* 

C0NM0N/C0H4/SIPPI 29.23) ,S I MPI 29.29) .Si *“*(29 .29) .SlPM(29.23),Pl.aMr 
l,SiPO(29.23).StMa(29.23).StOP(23.23>,SlOMI29.29)»SiOOI29.29) 
CQMM0N/CQM9/UCNTR. f£NTM. I WAULS .TWAULY »P ,(UNV .V2EPO* TSERO.EXPV.EXPU 
CGMNCN/COM*/ A* (29 .29 ) . ASJ ( 29 > . AS l ( 29) 

COMMON/COMT/ At PH AY (29) .BET AV( 29) .OEUTAr (29 1 .GAMMAY ( 29 1 . ALPHAS ( 29 ) . 

I 8 ETAS (29) .OELTASI 29) .GAMMAS (29) .V ISCK 29, 29 ). VI SC J (29, 29) 

CONMON/COM9/HC I GMT .-IDrH.HSLOPB, -SLOPE, OLOA. ASP* 

DIMENSION AU( 29. 29). 0(29.29) 

90 10 J-2.NJ 
00 10 1-2. NI 

10 *U(I ,4)*8Mail.J)«U(l»J) 

SUM-0 .0 
00 20 I«9 *NI 

20 SUM-SUH-SIPM( I-l .2 )/( 3.-8UI (-1 ,2) -PUI 1 .2) )/ GAMMAS ( I-l >-SIMM( 1 ,2 ) / 

1 (l.MU(I.2)*AUII-l.2) )/GAMNA2(I) 

02(2)-HHM(2.2)/GAMMA2(l)/8UI2»2)*l6.»SUM 
00 29 1*9 »Nl 
00 29 4-l.NJ 

29 G( I,J)-SIMN( l.J)/(2.*PU(l.4)*AU) I-l . J ) -AUI I , J-l ) )-SlPM( l- 1 . 0) / 

K2.-AUU-I.J) ♦AU( l.J )*«U(I-l, J-l) ) *SIPP( I— 1 * J— l )/(2.**U (1-1*4— 1)4* 
2Ull.J-l)4*U( l-l.J) )-SlMP( I * J-l )/( 2.-*U( l ,4- l )**U( I-l » ) ♦HU I 1 .4 1 ) 

00 19 4- l.NJ 
SUM-0 .0 

00 10 1-3, NI 
10 SUM-SUM-GI 1.9) 

SUM- SUM -S 1 4M( 2.4) /(l.**U(2 .4) *PU( 2, J-l ) )/GAMMAY(4 )*SIMP(2.J-l )/ 

1 (1.»*U(2»4-1)*«U12,J) )/GAMMAY(J-ll 

39 OZ(J)-SUN-l*. 

SUN-0. 

00 40 4-2. NJ 

40 SUM- SUM-02 (J) 

YN I-l ./(SUM-ASP* ) —0.9 
00 49 J-2.N4 
02(4)-0ZI JI-YN I 
49 2( 4)-2( J-l )*02(4) 

SUM-0. 

OC 47 4-3, N4 

47 SUH*SUM*$IMP(2,4-U/I3.-AU(2,4-1)-AU(2,4I)/GAMMAY( 4-l)-SlMM<2,4)/ 
l (3.*HU(2,J)*l»U(2,4-l) )/GAMMA2(J) 

OY ( 2) -SIMM (2,2 ) /GAMMAS (1)/*U(2»2)* 16. -SUM 

00 99 I-l, MI 

SUM-0. 

00 90 J-1.N4 
90 SUM-SUM-G(I.J) 
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UN [V 01 RELF 4 SE 2.0 20 ZV 0 V 04 TE « 91224 11 / 39/29 

$UM«SUM»SIM*( 1 . 2 )/l 3 .*RJ( I . 2 ) *RU< I-l» 2 >) /G 4 MM 4 Z ( 1 - 1 . 2 )/ 

l (3.MUI 1-1.2 ) ♦HU C 1.2) )/GAMNAZ( (-1) 

59 0 V ( 1 )*SUH* 16 . 

SUM. 0 . 

CO 60 l- 2 .Nl 

60 SUN-SUN-OV ( 1 1 
2 Nl-l./«SUW/ 4 SR*U— 0.5 
00 65 I- 2 .NI 
OVtU-OVC l)- 2 M 

65 ym«ru-u*oy(ii 

RETURN 

6 N 0 
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-vi (V ui RELEASE 2.0 V SC STY CATE > 11224 U/94/2E 

SUBROUTINE VSCSTYINP) 

CONNQN/COMl/’MI *9.29) .T( 25.29 1 .HI29.29) .PMO 1 29.291 »USHt 29.25) 
COHMON/CnM^r^ «*J1,NJ2.NJP.NI .Nil ,Nt2.M(».R«t2.N»CY,OK 
CONHGN/CON7/H?*«AVU9).EBTAYt29).OiLTAVI29».GANNAYt29).AlPMAj:t29> 
t BET A* (25) rO«UA2( 291 .GANNA2I2S) , VI SCI ( 29.291 .V tSC4(29.29) 

MMALL>M( 1.1) 

CC 44 ) J-UNJ 
OU 40 l-l.Nt 
V1SCIU. 4)0. 000019 
O VISC4I I (4)>0. 000019 
RETURN 
END 


* 
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sAN IV jI RELEASE 2.3 PflSOLV DATE • 8122* i/36/28 

SUBROUTINE PQSCLVI A l ,C l , it, AJ , CJ * 00, X , XX I 
CCHNON/CCHJ/NJ ,NJ 1 ,N J2 ,NJR , Nf » N 1 1 » N 1 2 » N l P ,N8L2 ,NBLY .OX 
0 IMENS ICMAII 25, 251 , CH 25.251 ,88(25.231 ,AJI 25 , 25 1 ,CJ 125, 25) 
l ,CD( 23,23),XU3,i:3),XX(23,25) .01 25) ,£( 25) ,f ( 23) .8(25) 

„C** ROUTINE SOLVES THE PENTAO tAGCNAL EQS. BV A-0- 1 PETHGO ••••••••••• 

00 50 t-l.NIP 
00 30 J-l.NJP 
50 XXI l,J)-X( I,,)) 

00 300 NN-2.MJ 

J 

00 130 1 ■ 2 , N I 

130 3( 1 1-301 I ,J)MX(l ,J-l)-CJ( l, < mXXII,J + l)»AJU ,J) 
£I2)-AI(2,J)/BB(2,J) 

PI 2) »0< 21/88 (2,3) 

00 UO 1-3, Nil 

DEN-BSl I , J )-C f II , J ) *v. U“l ) 

SI I) “All I,J)/OEN 

uo pin-iom^ci i i,j)*pi i-m/oEN 

XX(M(,J)-lO(Nn»P(N(ll*Ul (NI»J)*CI1NI,JI))/(88(NI,J)-G(N11)* 

1 (AUNI,J)KI(NI,J))) 

DO 120 N-I.NI2 
I-Nl-N 

120 XX( I , J) -El l)-XX(t*l«J)YP( 1) 

I-NN 

00 200 J-2.NJ 

200 01 4) «00< l , J ) ♦XX < I-l,J) Kl 1 1 ,J)tXX( l ♦ l . J ) *AI 1 I , J ) 

E ( 2) "AJ( 1 .21/381 (.2) 

P 1 2) -01 21/591 1,2) 

00 210 J-3.NJI 

0EN-6BI I,.J)-CJ< t,J)*€IJ-ll 

EUl-AOit ,J)/OEN 

213 P(J)-10IJ)>CJ( l,J)*P(J-l) )/0EN 

XXII,Na»-(0(NJ)+PINJl)*(AJlUNJ)M;JI l,NJ ) ) )/ IBBII ,NJ l-EINJ U- 
l t AJ 1 l,NJ)K.J< I.NJI)) 

OC 220 N-I.NJ2 
J-NJ-N 

223 XXII ,J)-S< J)*XXI 
300 CONTINUE 
HETUAN 
END 
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